package com.github.zyf.ioaj.algorithm.physics;

import com.github.zyf.ioaj.algorithm.AbstractAlgorithm;
import com.github.zyf.ioaj.annotation.internal.Algorithm;
import com.github.zyf.ioaj.annotation.internal.BasicParam;
import com.github.zyf.ioaj.annotation.internal.TempParam;
import com.github.zyf.ioaj.constant.FileType;
import com.github.zyf.ioaj.constant.ParamType;
import com.github.zyf.ioaj.export.*;
import com.github.zyf.ioaj.algorithm.Input;
import com.github.zyf.ioaj.algorithm.Output;

/**
 * <h1>Simulate Anneal Arithmetic</h1>
 * <p>模拟退火优化算法</p>
 * <p>默认求解最小值，如果需要求解最大值，将目标函数乘以-1即可，最后得到的值再次乘以-1就是得到的最优值（自变量取值不用乘以-1）</p>
 * <p>需要计算则调用{@link #go()}</p>
 * <p>需要跟踪每次迭代则调用{@link #tracker(String, String, String, FileType)}</p>
 * <p>需要生成报告则调用{@link #report(String, String, String, FileType)}</p>
 * <p>模拟退火算法基本思想如下：</p>
 * <p>
 * 假设有一温度很高的物体，其在自然降温过程中：
 * 当温度较高时，物体分子的热运动速率较大，分子较为活跃；
 * 当温度逐步下降时，物体分子的热运动速率较小，分子较为不活跃。
 * </p>
 * <p>
 * 我们将上述过程带入到一个函数求极值问题中，并假设自变量即为分子。
 * 那么，温度较高时自变量取值发生变化的幅度较大；温度较低时自变量取值发生变化的幅度较小。
 * </p>
 * <p>
 * 另外，当自变量取到一个较好的取值时，我们无条件接受取值；
 * 而当自变量取到一个不那么好甚至较差的取值时，我们以一定的概率接受取值。
 * 因为，取值较差的自变量也有可能经过变化变化到最优解附近。
 * 概率接受非较好取值，保证了算法跳出局部最优解避免陷入局部最优。
 * </p>
 * <p>
 * 具体计算方法如下：
 * <ul>
 *     <li>温度计算：Tc = alpha * Tc</li>
 *     <li>自变量计算：p = p + (r - 0.5) * Tc</li>
 *     <li>接受概率：采用Metropolis准则，metropolis_P = exp((p_old - p_new) / Tc)</li>
 * </ul>
 * 其中
 * <ul>
 *     <li>alpha：{@link #alpha}</li>
 *     <li>Tc：{@link #Tc}</li>
 *     <li>p：{@link #p}</li>
 *     <li>metropolis_P：{@link #metropolis_P}</li>
 *     <li>r：一个随机数</li>
 *     <li>p_old：上一次位置（旧值）</li>
 *     <li>p_new：当前位置（新值）</li>
 * </ul>
 *
 * @author zyf
 */
@Algorithm
public class SAA extends AbstractAlgorithm {

    /**
     * 分子个数
     */
    @BasicParam(type = ParamType.INT, desc = "分子个数")
    private int n;

    /**
     * 温度衰减系数，取值越接近1（不等于1），温度衰减越慢、求解精度越高，但求解时间越久。
     */
    @BasicParam(type = ParamType.DOUBLE, desc = "温度衰减系数")
    private double alpha;

    /**
     * 初始温度
     */
    @BasicParam(type = ParamType.DOUBLE, desc = "初始温度")
    private double T0;
    /**
     * 当前温度，将从{@link #T0}变化到{@link #Tf}
     */
    @TempParam(type = ParamType.DOUBLE, desc = "当前温度")
    private double Tc;
    /**
     * 结束温度
     */
    @BasicParam(type = ParamType.DOUBLE, desc = "结束温度")
    private double Tf;

    /**
     * 分子当前所在的位置
     */
    @TempParam(type = ParamType.DOUBLE_ARRAY_2, desc = "分子当前所在的位置")
    private double[][] p;

    /**
     * Metropolis准则计算出来的接受概率
     */
    @TempParam(type = ParamType.DOUBLE, desc = "TempParam")
    private double metropolis_P;

    /**
     * 对于当前的值是否接受（取决于Metropolis准则计算出来的接受概率 {@link #metropolis_P}）
     */
    @TempParam(type = ParamType.BOOL, desc = "对于当前的值是否接受")
    private boolean accept;

    /**
     * 当前全局最优值，double[input().varNum()]
     */
    @TempParam(type = ParamType.DOUBLE_ARRAY_1, desc = "当前全局最优值")
    private double[] x;

    /**
     * 自变量个数
     */
    private int D;
    /**
     * 自变量取值下限
     */
    private double[] lb;
    /**
     * 自变量取值下限
     */
    private double[] ub;

    /**
     * 抽象算法构造函数，需要提供输入量，即提前初始化好{@link Input}对象
     *
     * @param input {@link Input}
     */
    public SAA(Input input) {
        super(input);
        init();
    }

    /* 初始化部分 */

    public void init() {
        this.T0 = 100.0;
        this.Tc = T0;
        this.Tf = 0.01;
        this.alpha = 0.99;
        this.n = 50;

        this.D = getInput().getVarNum();
        this.lb = getInput().getLb();
        this.ub = getInput().getUb();

        this.p = new double[n][D];
        initP();
        initGbest();
    }

    public void initP() {
        int i, j;
        for (i = 0; i < n; i++) {
            for (j = 0; j < D; j++) {
                p[i][j] = lb[j] + Math.random() * (ub[j] - lb[j]);
            }
        }
    }

    public void initGbest() {
        int i;
        x = p[0].clone();
        for (i = 0; i < n; i++) {
            if (getInput().getFunction().apply(p[0]) < getInput().getFunction().apply(x)) x = p[i].clone();
        }
    }

    /* 计算部分 */
    /* 考虑到迭代次数和分子数目很多，数量级可能很大，不建议做形式上的优化 */

    /**
     * 每次迭代时扫描一次
     */
    public void scan() {
        int i, j;
        // 如果温度低于设定的终止温度，那么在逻辑上也就终止迭代了
        while (Tc > Tf) {
            // 扫描所有分子
            for (i = 0; i < n; i++) {
                double[] new_p = new double[D];
                // 扫描所有变量
                for (j = 0; j < D; j++) {
                    // 计算新解位置，随机数-0.5是为了保证正负分布均匀
                    new_p[j] = p[i][j] + (Math.random() - 0.5) * Tc;
                    // 防止越界
                    if (new_p[j] > ub[j]) new_p[j] = ub[j];
                    if (new_p[j] < lb[j]) new_p[j] = lb[j];
                }
                // 计算新值和旧值
                double new_val = getInput().getFunction().apply(new_p);
                double old_val = getInput().getFunction().apply(p[i]);
                // 更新群体最优解
                if (new_val < getInput().getFunction().apply(x)) x = new_p.clone();
                // 判断是否更新解
                metropolis(old_val, new_val);
                if (accept) p[i] = new_p.clone();
            }
            Tc *= alpha;
        }
    }

    /**
     * 生成最优解
     */
    private void generateOptimization() {
        // 获取最优值
        double optimalOutput = getInput().getFunction().apply(x);
        // 获取最优值时自变量取值
        double[] optimalInput = x.clone();
        setOutput(new Output(optimalOutput, optimalInput));
    }

    /**
     * 运行并计算结果
     */
    @Override
    public void go() {
        int t;
        for (t = 0; t < getIter(); t++) {
            scan();
            // 如果温度低于设定的终止温度，那么在逻辑上也就终止迭代了
            if (Tc < Tf) break;
        }
        generateOptimization();
    }

    /**
     * 运行并计算结果（按迭代逐步跟踪），如果目的是计算出结果，请使用{@link #go()}
     *
     * @param author   作者
     * @param message  备注
     * @param title    标题
     * @param fileType 文件类型 {@link FileType}
     * @return {@link Tracker}
     */
    @Override
    public Tracker tracker(String author, String message, String title, FileType fileType) {
        int t;
        Tracker tracker = new Tracker(this, author, message, title, fileType);
        tracker.startTracker();
        for (t = 0; t < getIter(); t++) {
            scan();
            tracker.trackOne(t + 1);
        }
        generateOptimization();
        tracker.endTracker();
        return tracker;
    }

    /**
     * 生成报告
     *
     * @param author   作者
     * @param message  备注
     * @param title    标题
     * @param fileType 文件类型 {@link FileType}
     * @return {@link Report}
     */
    @Override
    public Report report(String author, String message, String title, FileType fileType) {
        return new Report(this, author, message, title, fileType);
    }

    /**
     * Metropolis 准则
     *
     * @param old_val 旧值
     * @param new_val 新值
     */
    private void metropolis(double old_val, double new_val) {
        if (new_val < old_val) accept = true;
        else {
            // 计算概率
            metropolis_P = Math.exp((old_val - new_val) / Tc);
            // 根据概率判断是否接受
            accept = Math.random() < metropolis_P;
        }
    }

    /* Getter/Setter 部分 */

    public Integer getN() {
        return n;
    }
    public void setN(Integer n) {
        this.n = n;
    }
    public double getAlpha() {
        return alpha;
    }
    public void setAlpha(double alpha) {
        this.alpha = alpha;
    }
    public double getT0() {
        return T0;
    }
    public void setT0(double t0) {
        T0 = t0;
    }
    public double getTf() {
        return Tf;
    }
    public void setTf(double tf) {
        Tf = tf;
    }

    /* Getter部分 */

    public double getTc() {
        return Tc;
    }
    public double[][] getP() {
        return p;
    }
    public double getMetropolis_P() {
        return metropolis_P;
    }
    public boolean isAccept() {
        return accept;
    }
    public double[] getX() {
        return x;
    }
}
